#include <ctime>
#include <iostream>
#include <cmath>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <cstdlib>
#include <cstring>

#include <vector>
#include <string>
#include <random>

#include "tableaux_sections_efficaces.cpp"
#include "tableaux_energies.cpp"

using namespace std;

#define PI 3.14159265
#define electron_mass_c2 0.5109989461 * 1e6
#define N 1 
#define NB_E 100000
#define ENERGIE_INITIALE 100

struct data_res {
	int min;
	int max;
	int trouve;
	int niveauEnergie;
	double E;
	double sigma_sec_eff;
	double sigma_sec_eff_min;
	double sigma_sec_eff_max;
};



mt19937_64 gen(1729);
uniform_real_distribution<double> distribution(0.0, 1.0);



//Fonctions Communes ------------------------------------------------------------------

inline double frand(double min, double max);
inline double UniformRand();
inline double Interpolation(double Xa, double Ya, double Xb, double Yb, double X);
double QuadInterpolation(double e11, double e12, double e21, double e22, double xs11, double xs12, double xs21, double xs22, double t1, double t2, double t, double e);
inline int inside(double Wx, double Wy, double Wz, double x, double y, double z);


//Fonctions Collision Elastique --------------------------------------------------------
void getDirectionElastique(double E, double &Vx, double &Vy, double &Vz);
double getCosTheta(double E);
double getSigmaElastique(double E);
double GetLambdaElastique(double E);
double GetCosTheta(double E);
double GetPhi();
double gamma(double E);
double beta(double E);
double delta(double E);
double s(double E);

double DiffSigElastic(double E, double O);


//Fonctions Vibration ------------------------------------------------------------------
double getSigmaVibration(data_res &data);
int getIndiceEnergieVibration(data_res &data);
double getSigmaVibration(double E);
double GetLambdaVibration(double E);
int GetIndiceEnergieVibration(double E);


//Fonctions Excitation ------------------------------------------------------------------
double getSigmaExcitation(data_res &data);
int getIndiceEnergieExcitation(data_res &data);
double getSigmaExcitation(double E);
double GetLambdaExcitation(double E);
int GetIndiceEnergieExcitation(double E);


//Fonctions Ionisation ------------------------------------------------------------------
double getSigmaIonisation(data_res &data);
int getIndiceEnergieIonisation(data_res &data);
double getSigmaIonisation(double E);
double GetLambdaIonisation(double E);
int GetIndiceEnergieIonisation(double E);
void getDirectionElectronPrimaireSecondaire(const double &Ti, const double &Ts, double &Vx, double &Vy, double &Vz, double &Vxs, double &Vys, double &Vzs);
void GetDirectionElectronSecondaire(const double & Ti, const double & Ts, double & cosTheta, double & phi);
double QuadInterpolation(double e11, double e12, double e21, double e22, double xs11, double xs12, double xs21, double xs22, double t1, double t2, double t, double e);
double GetEnergieElectronSecondaire(double k, int indice);
double DifferentialCrossSection(double Ti, double Tt, int Indice);
void GetIndiceSigmaDiff(double Ti, int & min_inf, int & min_sup, int & max_inf, int & max_sup);


//Fonction Principale ---------------------------------------------------------------------------------------
int main()
{
	cout << "Simulation demarree ..." << endl;

	//Declaration des variables des fichiers de stockage de donnees ------------------------------------
	ofstream *outFile_xyz_attribut = new ofstream[N];
	ofstream *outFile_xyz_energie_elec_tous = new ofstream[N];

	stringstream sstm;

	int j = 0; //One simulation


	//Ouvrir les fichiers de stockage de donnees --------------------------------------------------------
	sstm.str("");
	sstm << "XYZ_Energies_All_" << j << ".dat";

	outFile_xyz_energie_elec_tous[j].open(sstm.str().c_str(), ios::binary);

	outFile_xyz_energie_elec_tous[j].width(2);
	outFile_xyz_energie_elec_tous[j] << setfill(' ') << setw(15) << showpoint << scientific;

	
	sstm.str("");
	sstm << "XYZ_attribut_" << j << ".dat";

	outFile_xyz_attribut[j].open(sstm.str().c_str(), ios::binary);

	outFile_xyz_attribut[j].width(2);
	outFile_xyz_attribut[j] << setw(15) << showpoint << scientific;
	//-----------------------------------------------------------------------------------------------

	//Declaration des variables de calcul
	int i, first, indice_energie, nb = NB_E, job_nb = N, k = 1;
	double Energie = ENERGIE_INITIALE;
	double  Xi, Yi, Zi, Xf, Yf, Zf, dist = 0, Vx, Vy, Vz, Ti, Td, Tf, Ts, Vxs, Vys, Vzs;//Xs = 0, Ys = 0, Zs = 0,
	double Wx = 1000000000, Wy = 1000000000, Wz = 1000000000;
	data_res dataVibration, dataExcitation, dataIonisation;


	double sigma_interaction[4], Lambda, somme_sigma, val_aleatoire;
	//lambda_interaction :   0 -> elastique, 1 -> vibration, 2 -> excitation, 3 -> ionisation
	int nb_interaction[4] = {0}, interaction, indice_sigma;

	double cutoff_energy[4] = { 8.22, 8.22, 8.22, 10.79 };
	//cutoff_energy :   0 -> simulation, 1 -> vibration, 2 -> excitation, 3 -> ionisation


	vector<double> TXs, TYs, TZs, TVxs, TVys, TVzs, Tts;

	//Declartion des variable de calcul de temps ---------------------------------------------------------
	clock_t start, end;
	double cpu_time_used;


	outFile_xyz_attribut[j] << "XYZ_Energies_All_" << j << ", NB Electron = " << NB_E << ", Energie initiale = " << ENERGIE_INITIALE << " eV, Wx = Wy = Wz = " << Wx << endl;

	//start timing ---------------------------------------------------------------------------------------
	start = clock();

	//Simulation pour NB_E electrons ---------------------------------------------------------------------
	for (i = 1; i <= nb; ++i)
	{
		//initialisation d'un nouveau electron ---------------------------------------
		Xi = Yi = Zi = 0;
		Xf = Yf = Zf = 0;
		Vx = 1;
		Vy = Vz = 0;
		Ti = Energie;
		first = 1;

		while (Ti >= cutoff_energy[0] && inside(Wx, Wy, Wz, Xf, Yf, Zf)) //Ti >= cutoff_energy && inside(Wx, Wy, Wz, Xf, Yf, Zf)
		{
			
			//Get Lambda et choix du type d'interaction -----------------------------------------------
			
			somme_sigma = 0;

			//calcul sigma Elastique
			
			somme_sigma = getSigmaElastique(Ti);
			sigma_interaction[0] = somme_sigma;
			
			//calcul sigma Vibration
			if (Ti >= cutoff_energy[1] && Ti <= 100.)
			{
				dataVibration.E = Ti;
				somme_sigma += getSigmaVibration(dataVibration);
			}
			sigma_interaction[1] = somme_sigma;

			//calcul sigma Excitation
			if (Ti >= cutoff_energy[2])
			{
				dataExcitation.E = Ti;
				somme_sigma += getSigmaExcitation(dataExcitation);
			}
			sigma_interaction[2] = somme_sigma;

			//calcul sigma Ionisation
			if (Ti >= cutoff_energy[3])
			{
				dataIonisation.E = Ti;
				somme_sigma += getSigmaIonisation(dataIonisation);
			}
			sigma_interaction[3] = somme_sigma;
			
			
			//Choix interaction
			val_aleatoire = frand(0., somme_sigma);

			for (indice_sigma = 0; indice_sigma < 4; ++indice_sigma)
				if (val_aleatoire < sigma_interaction[indice_sigma])
					break;

			interaction = indice_sigma;
			
			Lambda = -log(frand(0., 1.)) / somme_sigma; //libre parcour moyen = zetha / (1/lamda1 + 1/lambda2 + ...) 

			nb_interaction[interaction]++; //compter le nombre de chaque interaction


			switch (interaction)
			{
			case 0:  //Collision elastique ------------------------------------------------------------
					 //Calculer la nouvelle position de l'electron incident
				Xf = Lambda * Vx + Xi;
				Yf = Lambda * Vy + Yi;
				Zf = Lambda * Vz + Zi;

				getDirectionElastique(Ti, Vx, Vy, Vz);

				Xi = Xf;
				Yi = Yf;
				Zi = Zf;

				//A verifier---------------------------------------------------------------------------------------
				if(inside(Wx, Wy, Wz, Xf, Yf, Zf))
					outFile_xyz_energie_elec_tous[j] << Xf << "   " << Yf << "   " << Zf << "   " << Ti << "\n";
				//-------------------------------------------------------------------------------------------------
				break;

			case 1: //Vibration ---------------------------------------------------------------------------------
					//Get Energie deposee
				Td = energies_vibration[getIndiceEnergieVibration(dataVibration)];

				//Calculer la nouvelle position de l'electron incident
				Xf = Lambda * Vx + Xi;
				Yf = Lambda * Vy + Yi;
				Zf = Lambda * Vz + Zi;

				//Ajouter les nouvelles information au fichier de sortie
				outFile_xyz_energie_elec_tous[j] << Xf << "   " << Yf << "   " << Zf << "   " << Td << '\n';

				//Mise a jour des parametres de l'electron incident
				Xi = Xf;
				Yi = Yf;
				Zi = Zf;
				Ti = Ti - Td;
				break;

			case 2: //Excitation ------------------------------------------------------------------------------
					//Get Energie deposee
				Td = energies_excitation[getIndiceEnergieExcitation(dataExcitation)];

				//Il faut avancer l'electron
				Xf = Lambda * Vx + Xi;
				Yf = Lambda * Vy + Yi;
				Zf = Lambda * Vz + Zi;

				//Ajouter les nouvelles information au fichier de sortie
				outFile_xyz_energie_elec_tous[j] << Xf << "   " << Yf << "   " << Zf << "   " << Td << '\n';

				//Mise a jour des parametres de l'electron incident
				Xi = Xf;
				Yi = Yf;
				Zi = Zf;
				Ti = Ti - Td;
				break;

			case 3: //Ionisation -------------------------------------------------------------------------------
					// Get Energie deposee et get lambda electron incident
				indice_energie = getIndiceEnergieIonisation(dataIonisation);
				Td = energies_ionisation[indice_energie];


				//Calculer la nouvelle position de l'electron incident
				Xf = Lambda * Vx + Xi;
				Yf = Lambda * Vy + Yi;
				Zf = Lambda * Vz + Zi;

				//Ajouter les nouvelles information au fichier de sortie
				outFile_xyz_energie_elec_tous[j] << Xf << "   " << Yf << "   " << Zf << "   " << Td << '\n';

				//Get Energie de l'electron secondaire
				Ts = GetEnergieElectronSecondaire(Ti, indice_energie);

				getDirectionElectronPrimaireSecondaire(Ti, Ts, Vx, Vy, Vz, Vxs, Vys, Vzs);				


				//Energie restante de l'electron incident
				Tf = Ti - (Td + Ts);

				//Memoriser parametres electron secondaire pour simulation ulterieur :
				//Position Xf, Yf, Zf
				//Energie Ts
				//Direction Vxs, Vys, Vzs

				if (Ts >= cutoff_energy[0])
				{
					
					TXs.push_back(Xf);
					TYs.push_back(Yf);
					TZs.push_back(Zf);
					Tts.push_back(Ts);
					TVxs.push_back(Vxs);
					TVys.push_back(Vys);
					TVzs.push_back(Vzs);

				}
				

				//Mise a jour des parametres de l'electron incident
				Xi = Xf;
				Yi = Yf;
				Zi = Zf;
				Ti = Tf;

				break;
			}

			

			if (Ti < cutoff_energy[0])
			{
				//Avancer l'electron un dernier pas avec un collision elastique
				Lambda = GetLambdaElastique(Ti);

				Xf = Lambda * Vx + Xi;
				Yf = Lambda * Vy + Yi;
				Zf = Lambda * Vz + Zi;


				//Continuer la simulation avec un electron secondaire s'il y en a un
				if (!TXs.empty())//if (!TElecSec.empty())
				{
					
					Xi = TXs.back();
					Yi = TYs.back();
					Zi = TZs.back();
					Ti = Tts.back();
					Vx = TVxs.back();
					Vy = TVys.back();
					Vz = TVzs.back();
					first = 0;

					TXs.pop_back();
					TYs.pop_back();
					TZs.pop_back();
					Tts.pop_back();
					TVxs.pop_back();
					TVys.pop_back();
					TVzs.pop_back();

				}
			}

		}//end while (Ti >= cutoff_energy)



	}// end for (i = 1; i <= NB_E; i++)

	 //stop timing
	end = clock();

	cpu_time_used = ((double)(end - start)) / CLOCKS_PER_SEC;

	outFile_xyz_attribut[j] << endl << " CPU TIME USED : " << cpu_time_used << endl;

	cout << endl << " CPU TIME USED : " << cpu_time_used << endl;

	outFile_xyz_attribut[j] << "Nombre interaction : " << endl;
	outFile_xyz_attribut[j] << "             - Collision elastique : " << nb_interaction[0] << endl;
	outFile_xyz_attribut[j] << "             - Vibration : " << nb_interaction[1] << endl;
	outFile_xyz_attribut[j] << "             - Excitation : " << nb_interaction[2] << endl;
	outFile_xyz_attribut[j] << "             - Ionisation : " << nb_interaction[3] << endl;


	//Fermer les fichiers -------------------------------------------------------------------------------
	outFile_xyz_attribut[j].close();
	
	outFile_xyz_energie_elec_tous[j].close();


	return 0;
}

//Fonctions Communes ------------------------------------------------------------------
inline double Interpolation(double Xa, double Ya, double Xb, double Yb, double X)
{
	double Y, A, B;

	A = (Yb - Ya) / (Xb - Xa);
	B = Ya - A * Xa;

	Y = A * X + B;

	return Y;
}

inline double frand(double min, double max) {

	return distribution(gen) * (max - min) + min;
}

inline double UniformRand()
{
	return distribution(gen);
}

inline int inside(double Wx, double Wy, double Wz, double x, double y, double z)
{
	if (x >= -Wx && x <= Wx && y >= -Wy && y <= Wy && z >= -Wz && z <= Wz)
		return 1;
	else
		return 0;
}

//Fonctions Collision Elastique ------------------------------------------------------------------
void getDirectionElastique(double E, double &Vx, double &Vy, double &Vz)

{
	double Max;
	double cosTheta, Phi, y, diffSig;
	double U, V, W, NormV;

	if (E <= 200.)
	{
		double gamma, beta, delta;

		//Calcul de gamma --------------------------------------------------------------------
		double powE[6] = { 1, E, E*E, E*E*E, E*E*E*E, E*E*E*E*E };

		static double g[14] = { -1.7013, -1.48284, 0.6331, -0.10911, 8.353e-3, -2.388e-4, -3.32517, 0.10996, -4.5255e-3, 5.8372e-5, -2.4659e-7, 2.4775e-2, -2.96264e-5, -1.20655e-7 };

		if (E <= 10.)
			gamma = exp(g[0] + g[1] * powE[1] + g[2] * powE[2] + g[3] * powE[3] + g[4] * powE[4] + g[5] * powE[5]);
		else if (E <= 100.)
			gamma = exp(g[6] + g[7] * powE[1] + g[8] * powE[2] + g[9] * powE[3] + g[10] * powE[4]);
		else //if (E <= 200.)
			gamma = g[11] + g[12] * powE[1] + g[13] * powE[2];

		//Calcul de beta ----------------------------------------------------------------------
		static double b[5] = { 7.51525, -0.419122, 7.2017e-3, -4.646e-5, 1.02897e-7 };

		beta = exp(b[0] + b[1] * powE[1] + b[2] * powE[2] + b[3] * powE[3] + b[4] * powE[4]);

		//Calcul de delta ---------------------------------------------------------------------
		static double d[5] = { 2.9612, -0.26376, 4.307e-3, -2.6895e-5, 5.83505e-8 };
		delta = exp(d[0] + d[1] * powE[1] + d[2] * powE[2] + d[3] * powE[3] + d[4] * powE[4]);

		Max = 1. / (4.* gamma * gamma) + beta / ((2. + 2.* delta)*(2. + 2.* delta));

		do {
			cosTheta = distribution(gen) * 2. - 1.;
			y = distribution(gen) * Max;
			diffSig = ((1. / pow((1. + 2.*gamma - cosTheta), 2)) + (beta / pow((1. + 2.*delta + cosTheta), 2)));

		} while (y > diffSig);
	}
	else
	{
		double s;

		//Calcul de s -------------------------------------------------------------------------
		static double mc2 = 0.511e6;
		double sc = 1.64 - 0.0825*log(E);
		s = sc * 1.7e-5*pow(10., 2. / 3.)*mc2 / (E*(E / mc2 + 2.));

		Max = 1. / 4. / s / s;

		do {
			cosTheta = distribution(gen) * 2. - 1.;
			y = distribution(gen) * Max;
			diffSig = 1. / pow((1. + 2.*s - cosTheta), 2);;

		} while (y > diffSig);
	}


	Phi = 2. * PI * frand(0., 1.);

	//Calcul de la nouvelle direction

	if (1 - Vz * Vz <= 0) Vz = 0.9999; //Pour Eviter la division par 0

	U = cosTheta * Vx + (sqrt(1 - cosTheta * cosTheta) * (Vx * Vz * cos(Phi) - Vy * sin(Phi))) / sqrt(1 - Vz * Vz);
	V = cosTheta * Vy + (sqrt(1 - cosTheta * cosTheta) * (Vy * Vz * cos(Phi) + Vx * sin(Phi))) / sqrt(1 - Vz * Vz);
	W = cosTheta * Vz - (sqrt(1 - cosTheta * cosTheta) * sqrt(1 - Vz * Vz) * cos(Phi));

	NormV = sqrt(U * U + V * V + W * W);

	//Mise a jour des parameteres de l'electron incident
	Vx = U / NormV;
	Vy = V / NormV;
	Vz = W / NormV;
}

double getCosTheta(double E)
{
	double Max;
	double randcostheta, y, diffSig;

	if (E <= 200)
	{
		double gamma, beta, delta;

		//Calcul de gamma --------------------------------------------------------------------
		double powE[6] = { 1, E, E*E, E*E*E, E*E*E*E, E*E*E*E*E };

		static double g[14] = { -1.7013, -1.48284, 0.6331, -0.10911, 8.353e-3, -2.388e-4, -3.32517, 0.10996, -4.5255e-3, 5.8372e-5, -2.4659e-7, 2.4775e-2, -2.96264e-5, -1.20655e-7 };

		if (E <= 10.)
			gamma = exp(g[0] + g[1] * powE[1] + g[2] * powE[2] + g[3] * powE[3] + g[4] * powE[4] + g[5] * powE[5]);
		else if (E > 10. && E <= 100.)
			gamma = exp(g[6] + g[7] * powE[1] + g[8] * powE[2] + g[9] * powE[3] + g[10] * powE[4]);
		else if (E > 100. && E <= 200.)
			gamma = g[11] + g[12] * powE[1] + g[13] * powE[2];

		//Calcul de beta ----------------------------------------------------------------------
		static double b[5] = { 7.51525, -0.419122, 7.2017e-3, -4.646e-5, 1.02897e-7 };

		beta = exp(b[0] + b[1] * powE[1] + b[2] * powE[2] + b[3] * powE[3] + b[4] * powE[4]);

		//Calcul de delta ---------------------------------------------------------------------
		static double d[5] = { 2.9612, -0.26376, 4.307e-3, -2.6895e-5, 5.83505e-8 };
		delta = exp(d[0] + d[1] * powE[1] + d[2] * powE[2] + d[3] * powE[3] + d[4] * powE[4]);

		Max = 1. / (4.* gamma * gamma) + beta / ((2. + 2.* delta)*(2. + 2.* delta));

		do {
			randcostheta = frand(-1., 1.);
			y = frand(0., Max);
			
			diffSig = ((1. / pow((1. + 2.*gamma - randcostheta), 2)) + (beta / pow((1. + 2.*delta + randcostheta), 2)));

		} while (y > diffSig);
	}
	else
	{
		double s;

		//Calcul de s -------------------------------------------------------------------------
		static double mc2 = 0.511e6;
		double sc = 1.64 - 0.0825*log(E);
		s = sc * 1.7e-5*pow(10., 2. / 3.)*mc2 / (E*(E / mc2 + 2.));

		Max = 1. / 4. / s / s;

		do {
			randcostheta = frand(-1., 1.);
			y = frand(0., Max);
			
			diffSig = 1. / pow((1. + 2.*s - randcostheta), 2);;

		} while (y > diffSig);
	}

	return(randcostheta);
}

double getSigmaElastique(double E)
{
	static double mc2 = 0.511e6;
	static double Z = 10.;
	static double e = 1.;
	static double fourpiEps = 4.*3.1415*5.52635e-2;//eps is in q2/eV.nm
	double length = (e*e * (E + mc2)) / (fourpiEps * E * (E + 2. * mc2)); //fourpiEps*1e9 in q2/eV.m
	double R = Z * (Z + 1) * length * length;

	double sc = 1.64 - 0.0825*log(E);
	double Term = sc * 1.7e-5*pow(10., 2. / 3.)*mc2 / (E*(E / mc2 + 2.));
	double sigma = PI * R / Term / (Term + 1.);

	return (sigma * 33.43);//nanometer
}

double GetLambdaElastique(double E)
{
	double mc2 = 0.511e6;
	double Z = 10.;
	double e = 1;
	double fourpiEps = 4.*3.1415*5.52635e-2;//eps is in q2/eV.nm
	double length = (e*e * (E + mc2)) / (fourpiEps * E * (E + 2. * mc2)); //fourpiEps*1e9 in q2/eV.m
	double R = Z * (Z + 1) * length * length;

	double sc = 1.64 - 0.0825*log(E);
	double Term = sc * 1.7e-5*pow(10., 2. / 3.)*mc2 / (E*(E / mc2 + 2.));
	double sigma = PI * R / Term / (Term + 1.);
	double calculated_MFP = 1. / (sigma*3.343 * 10);//nanometer
	

	return (sigma*3.343 * 10);
}

double DiffSigElastic(double E, double O)
{
	double S;
	

	if (E <= 200.)S = ((1. / pow((1. + 2.*gamma(E) - O), 2)) + (beta(E) / pow((1. + 2.*delta(E) + O), 2)));
	else if (E>200.)S = 1. / pow((1. + 2.*s(E) - O), 2);

	return(S);
}

double GetCosTheta(double E)
{

	double Max;
	if (E <= 200)Max = 1. / (4.*gamma(E)*gamma(E)) + beta(E) / ((2. + 2.*delta(E))*(2. + 2.*delta(E)));
	else Max = 1. / 4. / s(E) / s(E);

	double randcostheta;
	double y;

	do {
		randcostheta = frand(-1., 1.);
		y = frand(0., Max);
		

	} while (y>DiffSigElastic(E, randcostheta));



	return(randcostheta);


}

double GetPhi()
{
	return frand(0., 2. * PI);
}

double s(double E)
{
	double mc2 = 0.511e6;
	double sc = 1.64 - 0.0825*log(E);
	double Term = sc * 1.7e-5*pow(10., 2. / 3.)*mc2 / (E*(E / mc2 + 2.));
	return(Term);
}

double gamma(double E)
{
	double g[14] = { -1.7013, -1.48284, 0.6331, -0.10911, 8.353e-3, -2.388e-4, -3.32517, 0.10996, -4.5255e-3, 5.8372e-5, -2.4659e-7, 2.4775e-2, -2.96264e-5, -1.20655e-7 };
	if (E <= 10.)return(exp(g[0] + g[1] * E + g[2] * E*E + g[3] * E*E*E + g[4] * pow(E, 4) + g[5] * pow(E, 5)));
	else if (E > 10. && E <= 100.)return(exp(g[6] + g[7] * E + g[8] * E*E + g[9] * E*E*E + g[10] * pow(E, 4)));
	else if (E > 100. && E <= 200.)return((g[11] + g[12] * E + g[13] * E*E));
	else
	{
		cout << "Erreur !!! gamma utiliser pour une valeur > 200 eV" << endl;
		exit(1);
	}

}

double delta(double E)
{
	double d[5] = { 2.9612, -0.26376, 4.307e-3, -2.6895e-5, 5.83505e-8 };
	return(exp(d[0] + d[1] * E + d[2] * E*E + d[3] * pow(E, 3) + d[4] * pow(E, 4)));
}

double beta(double E)
{
	double b[5] = { 7.51525, -0.419122, 7.2017e-3, -4.646e-5, 1.02897e-7 };
	return(exp(b[0] + b[1] * E + b[2] * E*E + b[3] * pow(E, 3) + b[4] * pow(E, 4)));
}

//Fonctions Vibration -------------------------------------------------------------------
double getSigmaVibration(data_res &data)
{
	int min, max, trouve;
	double sigma;

	min = 0;
	max = 37;
	trouve = -1;

	while (min <= max && trouve == -1)
	{
		if (e_Vib_Sanche[(max + min) / 2][0] == data.E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (e_Vib_Sanche[(max + min) / 2][0]> data.E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;

	}

	data.min = min;
	data.max = max;
	data.trouve = trouve;

	if (trouve >= 0)
		sigma = 1. / e_Vib_Sanche[trouve][11];
	else
		sigma = 1. / Interpolation(e_Vib_Sanche[max][0], e_Vib_Sanche[max][11], e_Vib_Sanche[min][0], e_Vib_Sanche[min][11], data.E);

	return sigma;
}

int getIndiceEnergieVibration(data_res &data)
{
	int i;
	double somme, sigma_rand;


	if (data.trouve >= 0)
	{
		sigma_rand = frand(0, e_Vib_Sanche[data.trouve][10]);
		for (i = 1, somme = 0; i <= 9; ++i)
		{
			somme += e_Vib_Sanche[data.trouve][i];

			if (somme >= sigma_rand || i == 9)
				break;
		}

	}
	else
	{
		sigma_rand = frand(0, Interpolation(e_Vib_Sanche[data.max][0], e_Vib_Sanche[data.max][10], e_Vib_Sanche[data.min][0], e_Vib_Sanche[data.min][10], data.E));
		for (i = 1, somme = 0; i <= 9; ++i)
		{
			somme += Interpolation(e_Vib_Sanche[data.max][0], e_Vib_Sanche[data.max][i], e_Vib_Sanche[data.min][0], e_Vib_Sanche[data.min][i], data.E);

			if (somme >= sigma_rand || i == 9)
				break;
		}

	}


	return i - 1;
}

double getSigmaVibration(double E)
{
	int min, max, trouve;
	double sigma;

	min = 0;
	max = 37;
	trouve = -1;



	while (min <= max && trouve == -1)
	{
		if (e_Vib_Sanche[(max + min) / 2][0] == E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (e_Vib_Sanche[(max + min) / 2][0]>E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;

	}

	if (trouve >= 0)
		sigma = 1. / e_Vib_Sanche[trouve][11];
	else
		sigma = 1. / Interpolation(e_Vib_Sanche[max][0], e_Vib_Sanche[max][11], e_Vib_Sanche[min][0], e_Vib_Sanche[min][11], E);

	return sigma;
}

double GetLambdaVibration(double E)
{
	int min, max, trouve;
	double lambda;

	min = 0;
	max = 37;
	trouve = -1;



	while (min <= max && trouve == -1)
	{
		if (e_Vib_Sanche[(max + min) / 2][0] == E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (e_Vib_Sanche[(max + min) / 2][0]>E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;

	}


	if (trouve >= 0)
		lambda = e_Vib_Sanche[trouve][11];
	else
		lambda = Interpolation(e_Vib_Sanche[max][0], e_Vib_Sanche[max][11], e_Vib_Sanche[min][0], e_Vib_Sanche[min][11], E);

	return lambda;
}

int GetIndiceEnergieVibration(double E)
{
	int min, max, i, trouve;
	double somme, sigma_rand;

	min = 0;
	max = 37;
	trouve = -1;

	while (min <= max && trouve == -1)
	{
		if (e_Vib_Sanche[(max + min) / 2][0] == E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (e_Vib_Sanche[(max + min) / 2][0]>E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;
	}

	if (trouve >= 0)
	{
		sigma_rand = frand(0, e_Vib_Sanche[trouve][10]);
		for (i = 1, somme = 0; i <= 9; i++)
		{
			somme += e_Vib_Sanche[trouve][i];

			if (somme >= sigma_rand || i == 9)
				break;
		}

	}
	else
	{
		sigma_rand = frand(0, Interpolation(e_Vib_Sanche[max][0], e_Vib_Sanche[max][10], e_Vib_Sanche[min][0], e_Vib_Sanche[min][10], E));
		for (i = 1, somme = 0; i <= 9; i++)
		{
			somme += Interpolation(e_Vib_Sanche[max][0], e_Vib_Sanche[max][i], e_Vib_Sanche[min][0], e_Vib_Sanche[min][i], E);

			if (somme >= sigma_rand || i == 9)
				break;
		}

	}


	return i - 1;

}

//Fonctions Excitation ------------------------------------------------------------------
double getSigmaExcitation(data_res &data)
{
	int min, max, trouve, niveauEnergie;
	double sigma_sec_eff, sigma_sec_eff_min, sigma_sec_eff_max; 

	min = 0;
	max = 82;
	trouve = -1;



	while (min <= max && trouve == -1)
	{
		if (sigma_excitation_e_born[(max + min) / 2][0] == data.E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (sigma_excitation_e_born[(max + min) / 2][0] > data.E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;
	}

	data.min = min;
	data.max = max;
	data.trouve = trouve;

	if (trouve >= 0)
	{
		if (data.E >= energies_excitation[4])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3] + sigma_excitation_e_born[trouve][4] + sigma_excitation_e_born[trouve][5];
			niveauEnergie = 4;
		}
		else if (data.E >= energies_excitation[3])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3] + sigma_excitation_e_born[trouve][4];
			niveauEnergie = 3;
		}
		else if (data.E >= energies_excitation[2])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3];
			niveauEnergie = 2;
		}
		else if (data.E >= energies_excitation[1])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2];
			niveauEnergie = 1;
		}
		else if (data.E >= energies_excitation[0])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1];
			niveauEnergie = 0;
		}
		else
		{
			cout << "Erreur getSigmaExcitation " << endl;
			exit(1);
		}
		
		data.niveauEnergie = niveauEnergie;
		data.sigma_sec_eff = sigma_sec_eff;
	}
	else
	{
		if (data.E >= energies_excitation[4])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3] + sigma_excitation_e_born[min][4] + sigma_excitation_e_born[min][5];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3] + sigma_excitation_e_born[max][4] + sigma_excitation_e_born[max][5];
			niveauEnergie = 4;
		}
		else if (data.E >= energies_excitation[3])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3] + sigma_excitation_e_born[min][4];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3] + sigma_excitation_e_born[max][4];
			niveauEnergie = 3;
		}
		else if (data.E >= energies_excitation[2])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3];
			niveauEnergie = 2;
		}
		else if (data.E >= energies_excitation[1])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2];
			niveauEnergie = 1;
		}
		else if (data.E >= energies_excitation[0])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1];
			niveauEnergie = 0;
		}
		else
		{
			cout << "Erreur getSigmaExcitation" << endl;
			exit(1);
		}
		sigma_sec_eff = Interpolation(sigma_excitation_e_born[max][0], sigma_sec_eff_max, sigma_excitation_e_born[min][0], sigma_sec_eff_min, data.E);
		
		data.niveauEnergie = niveauEnergie;
		data.sigma_sec_eff_min = sigma_sec_eff_min;
		data.sigma_sec_eff_max = sigma_sec_eff_max;
	}

	return sigma_sec_eff / 1000;
}

int getIndiceEnergieExcitation(data_res &data)
{
	int i;
	double somme, sigma_rand;


	if (data.trouve >= 0)
	{
		if (data.niveauEnergie == 4)
			sigma_rand = frand(0, sigma_excitation_e_born[data.trouve][6]);
		else
			sigma_rand = frand(0, data.sigma_sec_eff);

		for (i = 1, somme = 0; i <= data.niveauEnergie; ++i)
		{
			somme += sigma_excitation_e_born[data.trouve][i];

			if (somme >= sigma_rand)
				break;
		}

	}
	else
	{
		if (data.niveauEnergie == 4)
			sigma_rand = frand(0, Interpolation(sigma_excitation_e_born[data.max][0], sigma_excitation_e_born[data.max][6], sigma_excitation_e_born[data.min][0], sigma_excitation_e_born[data.min][6], data.E));
		else
			sigma_rand = frand(0, Interpolation(sigma_excitation_e_born[data.max][0], data.sigma_sec_eff_max, sigma_excitation_e_born[data.min][0], data.sigma_sec_eff_min, data.E));

		for (i = 1, somme = 0; i <= data.niveauEnergie; ++i)
		{
			somme += Interpolation(sigma_excitation_e_born[data.max][0], sigma_excitation_e_born[data.max][i], sigma_excitation_e_born[data.min][0], sigma_excitation_e_born[data.min][i], data.E);

			if (somme >= sigma_rand)
				break;
		}

	}

	return i - 1;
}

double getSigmaExcitation(double E)
{
	int min, max, trouve;
	double sigma_sec_eff, sigma_sec_eff_min, sigma_sec_eff_max; 

	min = 0;
	max = 82;
	trouve = -1;



	while (min <= max && trouve == -1)
	{
		if (sigma_excitation_e_born[(max + min) / 2][0] == E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (sigma_excitation_e_born[(max + min) / 2][0]>E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;
	}


	if (trouve >= 0)
	{
		if (E >= energies_excitation[4])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3] + sigma_excitation_e_born[trouve][4] + sigma_excitation_e_born[trouve][5];
		}
		else if (E >= energies_excitation[3])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3] + sigma_excitation_e_born[trouve][4];
		}
		else if (E >= energies_excitation[2])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3];
		}
		else if (E >= energies_excitation[1])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2];
		}
		else if (E >= energies_excitation[0])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1];
		}
		else
		{
			cout << "Erreur getSigmaExcitation " << endl;
			exit(1);
		}
	}
	else
	{
		if (E >= energies_excitation[4])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3] + sigma_excitation_e_born[min][4] + sigma_excitation_e_born[min][5];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3] + sigma_excitation_e_born[max][4] + sigma_excitation_e_born[max][5];
		}
		else if (E >= energies_excitation[3])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3] + sigma_excitation_e_born[min][4];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3] + sigma_excitation_e_born[max][4];
		}
		else if (E >= energies_excitation[2])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3];
		}
		else if (E >= energies_excitation[1])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2];
		}
		else if (E >= energies_excitation[0])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1];
		}
		else
		{
			cout << "Erreur getSigmaExcitation " << endl;
			exit(1);
		}
		sigma_sec_eff = Interpolation(sigma_excitation_e_born[max][0], sigma_sec_eff_max, sigma_excitation_e_born[min][0], sigma_sec_eff_min, E);
	}

	return sigma_sec_eff / 1000;
}

double GetLambdaExcitation(double E)
{
	int min, max, trouve;
	double sigma_sec_eff, sigma_sec_eff_min, sigma_sec_eff_max; //lambda, 

	min = 0;
	max = 82;
	trouve = -1;



	while (min <= max && trouve == -1)
	{
		if (sigma_excitation_e_born[(max + min) / 2][0] == E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (sigma_excitation_e_born[(max + min) / 2][0]>E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;

	}


	if (trouve >= 0)
	{
		if (E >= energies_excitation[4])
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3] + sigma_excitation_e_born[trouve][4] + sigma_excitation_e_born[trouve][5];
		else if (E >= energies_excitation[3])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3] + sigma_excitation_e_born[trouve][4];
		}
		else if (E >= energies_excitation[2])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3];
		}
		else if (E >= energies_excitation[1])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2];
		}
		else if (E >= energies_excitation[0])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1];
		}
		else
		{
			cout << "Erreur GetLambdaExcitation  " << endl;
			exit(1);
		}
	}
	else
	{
		if (E >= energies_excitation[4])
		{
			
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3] + sigma_excitation_e_born[min][4] + sigma_excitation_e_born[min][5];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3] + sigma_excitation_e_born[max][4] + sigma_excitation_e_born[max][5];
			
		}
		else if (E >= energies_excitation[3])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3] + sigma_excitation_e_born[min][4];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3] + sigma_excitation_e_born[max][4];

		}
		else if (E >= energies_excitation[2])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3];

		}
		else if (E >= energies_excitation[1])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2];
		}
		else if (E >= energies_excitation[0])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1];
		}
		else
		{
			cout << "Erreur GetLambdaExcitation" << endl;
			exit(1);
		}
		sigma_sec_eff = Interpolation(sigma_excitation_e_born[max][0], sigma_sec_eff_max, sigma_excitation_e_born[min][0], sigma_sec_eff_min, E);
	}

	return sigma_sec_eff / 1000;
}

int GetIndiceEnergieExcitation(double E)
{
	int min, max, i, trouve, niveauEnergie;
	double somme, sigma_rand, sigma_sec_eff, sigma_sec_eff_min, sigma_sec_eff_max;

	min = 0;
	max = 82;
	trouve = -1;

	while (min <= max && trouve == -1)
	{
		if (sigma_excitation_e_born[(max + min) / 2][0] == E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (sigma_excitation_e_born[(max + min) / 2][0]>E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;
	}

	if (trouve >= 0)
	{
		if (E >= energies_excitation[4])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3] + sigma_excitation_e_born[trouve][4] + sigma_excitation_e_born[trouve][5];
			niveauEnergie = 4;
		}
		else if (E >= energies_excitation[3])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3] + sigma_excitation_e_born[trouve][4];
			niveauEnergie = 3;
		}
		else if (E >= energies_excitation[2])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2] + sigma_excitation_e_born[trouve][3];
			niveauEnergie = 2;
		}
		else if (E >= energies_excitation[1])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1] + sigma_excitation_e_born[trouve][2];
			niveauEnergie = 1;
		}
		else if (E >= energies_excitation[0])
		{
			sigma_sec_eff = sigma_excitation_e_born[trouve][1];
			niveauEnergie = 0;
		}
		else
		{
			cout << "Erreur GetIndiceEnergieExcitation" << endl;
			exit(1);
		}

		sigma_rand = frand(0, sigma_sec_eff);
		for (i = 1, somme = 0; i <= niveauEnergie; i++)
		{
			somme += sigma_excitation_e_born[trouve][i];

			if (somme >= sigma_rand)
				break;
		}

	}
	else
	{
		if (E >= energies_excitation[4])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3] + sigma_excitation_e_born[min][4] + sigma_excitation_e_born[min][5];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3] + sigma_excitation_e_born[max][4] + sigma_excitation_e_born[max][5];

			niveauEnergie = 4;
		}
		else if (E >= energies_excitation[3])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3] + sigma_excitation_e_born[min][4];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3] + sigma_excitation_e_born[max][4];
			niveauEnergie = 3;
		}
		else if (E >= energies_excitation[2])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2] + sigma_excitation_e_born[min][3];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2] + sigma_excitation_e_born[max][3];
			niveauEnergie = 2;
		}
		else if (E >= energies_excitation[1])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1] + sigma_excitation_e_born[min][2];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1] + sigma_excitation_e_born[max][2];
			niveauEnergie = 1;
		}
		else if (E >= energies_excitation[0])
		{
			sigma_sec_eff_min = sigma_excitation_e_born[min][1];
			sigma_sec_eff_max = sigma_excitation_e_born[max][1];
			niveauEnergie = 0;
		}
		else
		{
			cout << "Erreur GetIndiceEnergieExcitation " << endl;
			exit(1);
		}

		sigma_rand = frand(0, Interpolation(sigma_excitation_e_born[max][0], sigma_sec_eff_max, sigma_excitation_e_born[min][0], sigma_sec_eff_min, E));
		for (i = 1, somme = 0; i <= niveauEnergie; i++)
		{
			somme += Interpolation(sigma_excitation_e_born[max][0], sigma_excitation_e_born[max][i], sigma_excitation_e_born[min][0], sigma_excitation_e_born[min][i], E);

			if (somme >= sigma_rand)
				break;
		}

	}

	return i - 1;
}



//Fonctions Ionisation ------------------------------------------------------------------
double getSigmaIonisation(data_res &data)
{
	int min, max, trouve, niveauEnergie;
	double sigma_sec_eff = 0, sigma_sec_eff_min = 0, sigma_sec_eff_max = 0;

	min = 0;
	max = 82;
	trouve = -1;

	while (min <= max && trouve == -1)
	{
		if (sigma_ionisation_e_born[(max + min) / 2][0] == data.E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (sigma_ionisation_e_born[(max + min) / 2][0] > data.E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;

	}

	data.min = min;
	data.max = max;
	data.trouve = trouve;

	if (trouve >= 0)
	{
		if (data.E >= energies_ionisation[4])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3] + sigma_ionisation_e_born[trouve][4] + sigma_ionisation_e_born[trouve][5];
			niveauEnergie = 4;
		}
		else if (data.E >= energies_ionisation[3])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3] + sigma_ionisation_e_born[trouve][4];
			niveauEnergie = 3;
		}
		else if (data.E >= energies_ionisation[2])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3];
			niveauEnergie = 2;
		}
		else if (data.E >= energies_ionisation[1])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2];
			niveauEnergie = 1;
		}
		else if (data.E >= energies_ionisation[0])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1];
			niveauEnergie = 0;
		}
		else
		{
			cout << "Erreur : Section efficace = 0  getSigmaIonisation  " << endl;
			exit(1);
		}
		data.niveauEnergie = niveauEnergie;
		data.sigma_sec_eff = sigma_sec_eff;

	}
	else
	{
		if (data.E >= energies_ionisation[4])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3] + sigma_ionisation_e_born[min][4] + sigma_ionisation_e_born[min][5];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3] + sigma_ionisation_e_born[max][4] + sigma_ionisation_e_born[max][5];
			niveauEnergie = 4;
		}
		else if (data.E >= energies_ionisation[3])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3] + sigma_ionisation_e_born[min][4];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3] + sigma_ionisation_e_born[max][4];
			niveauEnergie = 3;
		}
		else if (data.E >= energies_ionisation[2])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3];
			niveauEnergie = 2;
		}
		else if (data.E >= energies_ionisation[1])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2];
			niveauEnergie = 1;
		}
		else if (data.E >= energies_ionisation[0])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1];
			niveauEnergie = 0;
		}
		else
		{
			cout << " Erreur : section efficace = 0  GetLambdaIonisation " << endl;
			exit(1);
		}

		sigma_sec_eff = Interpolation(sigma_ionisation_e_born[max][0], sigma_sec_eff_max, sigma_ionisation_e_born[min][0], sigma_sec_eff_min, data.E);

		data.niveauEnergie = niveauEnergie;
		data.sigma_sec_eff_min = sigma_sec_eff_min;
		data.sigma_sec_eff_max = sigma_sec_eff_max;
	}

	return sigma_sec_eff / 1000;
}

int getIndiceEnergieIonisation(data_res &data)
{
	int i;
	double somme, sigma_rand;

	

	if (data.trouve >= 0)
	{
		if(data.niveauEnergie == 4)
			sigma_rand = frand(0, sigma_ionisation_e_born[data.trouve][6]);
		else
			sigma_rand = frand(0, data.sigma_sec_eff);
		
		for (i = 1, somme = 0; i <= data.niveauEnergie; ++i)
		{
			somme += sigma_ionisation_e_born[data.trouve][i];

			if (somme >= sigma_rand)
				break;
		}
	}
	else
	{
		if (data.niveauEnergie == 4)
			sigma_rand = frand(0, Interpolation(sigma_ionisation_e_born[data.max][0], sigma_ionisation_e_born[data.max][6], sigma_ionisation_e_born[data.min][0], sigma_ionisation_e_born[data.min][6], data.E));
		else
			sigma_rand = frand(0, Interpolation(sigma_ionisation_e_born[data.max][0], data.sigma_sec_eff_max, sigma_ionisation_e_born[data.min][0], data.sigma_sec_eff_min, data.E));
		
		for (i = 1, somme = 0; i <= data.niveauEnergie; ++i)
		{
			somme += Interpolation(sigma_ionisation_e_born[data.max][0], sigma_ionisation_e_born[data.max][i], sigma_ionisation_e_born[data.min][0], sigma_ionisation_e_born[data.min][i], data.E);

			if (somme >= sigma_rand)
				break;
		}
	}


	return (i - 1);

}

double getSigmaIonisation(double E)
{
	int min, max, trouve;
	double sigma_sec_eff = 0, sigma_sec_eff_min = 0, sigma_sec_eff_max = 0;

	min = 0;
	max = 82;
	trouve = -1;

	while (min <= max && trouve == -1)
	{
		if (sigma_ionisation_e_born[(max + min) / 2][0] == E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (sigma_ionisation_e_born[(max + min) / 2][0]>E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;

	}

	if (trouve >= 0)
	{
		if (E >= energies_ionisation[4])
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3] + sigma_ionisation_e_born[trouve][4] + sigma_ionisation_e_born[trouve][5];
		else if (E >= energies_ionisation[3])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3] + sigma_ionisation_e_born[trouve][4];
		}
		else if (E >= energies_ionisation[2])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3];
		}
		else if (E >= energies_ionisation[1])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2];
		}
		else if (E >= energies_ionisation[0])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1];
		}
		else
		{
			cout << "Erreur : Section efficace = 0" << endl;
			exit(1);
		}

	}
	else
	{
		if (E >= energies_ionisation[4])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3] + sigma_ionisation_e_born[min][4] + sigma_ionisation_e_born[min][5];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3] + sigma_ionisation_e_born[max][4] + sigma_ionisation_e_born[min][5];
		}
		else if (E >= energies_ionisation[3])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3] + sigma_ionisation_e_born[min][4];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3] + sigma_ionisation_e_born[max][4];
		}
		else if (E >= energies_ionisation[2])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3];
		}
		else if (E >= energies_ionisation[1])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2];
		}
		else if (E >= energies_ionisation[0])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1];
		}
		else
		{
			cout << " Erreur : section efficace = 0 " << endl;
			exit(1);
		}

		sigma_sec_eff = Interpolation(sigma_ionisation_e_born[max][0], sigma_sec_eff_max, sigma_ionisation_e_born[min][0], sigma_sec_eff_min, E);
	}

	return sigma_sec_eff / 1000;
}

void getDirectionElectronPrimaireSecondaire(const double &Ti, const double &Ts, double &Vx, double &Vy, double &Vz, double &Vxs, double &Vys, double &Vzs)
{
	double PhiSecondaire, cosThetaSecondaire;
	double U, V, W;

	PhiSecondaire = 2. * PI * UniformRand();

	if (Ts < 50.)
		cosThetaSecondaire = (2.*UniformRand()) - 1.;
	else if (Ts <= 200.)
	{
		if (UniformRand() <= 0.1)
			cosThetaSecondaire = (2.*UniformRand()) - 1.;
		else
			cosThetaSecondaire = UniformRand()*(sqrt(2.) / 2.);
	}
	else
	{
		double sin2O = (1. - Ts / Ti) / (1. + Ts / (2.*electron_mass_c2));
		cosThetaSecondaire = sqrt(1. - sin2O);
	}


	U = sqrt(1 - cosThetaSecondaire * cosThetaSecondaire) * cos(PhiSecondaire);
	V = sqrt(1 - cosThetaSecondaire * cosThetaSecondaire) * sin(PhiSecondaire);
	W = cosThetaSecondaire;


	double up = Vx * Vx + Vy * Vy;

	if (up)
	{
		up = sqrt(up);
		double px = U, py = V, pz = W;

		U = (Vx * Vz * px - Vy * py + Vx * up * pz) / up;
		V = (Vy * Vz * px + Vx * py + Vy * up * pz) / up;
		W = (Vz * Vz * px - px + Vz * up * pz) / up;
	}
	else if (Vz < 0)
	{
		U = -U;
		W = -W;
	}


	Vxs = U;
	Vys = V;
	Vzs = W;

	double Psecondaire = sqrt(Ts * (Ts + 2 * electron_mass_c2));
	double TotalEnergy = Ti + electron_mass_c2;
	double Psquare = Ti * (TotalEnergy + electron_mass_c2);
	double TotalMomentum = sqrt(Psquare);

	double Px = TotalMomentum * Vx - Psecondaire * U;
	double Py = TotalMomentum * Vy - Psecondaire * V;
	double Pz = TotalMomentum * Vz - Psecondaire * W;

	double Pfinal = sqrt(Px * Px + Py * Py + Pz * Pz);
	Vx = Px / Pfinal;
	Vy = Py / Pfinal;
	Vz = Pz / Pfinal;

}

double GetLambdaIonisation(double E)
{
	int min, max, trouve;
	double sigma_sec_eff = 0, sigma_sec_eff_min = 0, sigma_sec_eff_max = 0; //lambda,

	min = 0;
	max = 82;
	trouve = -1;

	while (min <= max && trouve == -1)
	{
		if (sigma_ionisation_e_born[(max + min) / 2][0] == E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (sigma_ionisation_e_born[(max + min) / 2][0]>E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;

	}

	if (trouve >= 0)
	{
		if (E >= energies_ionisation[4])
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3] + sigma_ionisation_e_born[trouve][4] + sigma_ionisation_e_born[trouve][5];
		else if (E >= energies_ionisation[3])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3] + sigma_ionisation_e_born[trouve][4];
		}
		else if (E >= energies_ionisation[2])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3];

		}
		else if (E >= energies_ionisation[1])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2];

		}
		else if (E >= energies_ionisation[0])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1];

		}
		else
		{

			cout << "Erreur : Section efficace = 0 " << endl;
			exit(1);
		}
	}
	else
	{
		if (E >= energies_ionisation[4])
		{
			
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3] + sigma_ionisation_e_born[min][4] + sigma_ionisation_e_born[min][5];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3] + sigma_ionisation_e_born[max][4] + sigma_ionisation_e_born[min][5];
		}
		else if (E >= energies_ionisation[3])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3] + sigma_ionisation_e_born[min][4];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3] + sigma_ionisation_e_born[max][4];
			
		}
		else if (E >= energies_ionisation[2])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3];
			
		}
		else if (E >= energies_ionisation[1])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2];
			
		}
		else if (E >= energies_ionisation[0])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1];
			
		}
		else
		{
			cout << " Erreur : section efficace = 0 GetLambdaIonisation  " << endl;
			exit(1);
		}

		
		sigma_sec_eff = Interpolation(sigma_ionisation_e_born[max][0], sigma_sec_eff_max, sigma_ionisation_e_born[min][0], sigma_sec_eff_min, E);
	}


	return sigma_sec_eff / 1000;
}

int GetIndiceEnergieIonisation(double E)
{
	int min, max, i, trouve, temp_indice = -1;
	double somme, sigma_rand, sigma_sec_eff, sigma_sec_eff_min = 0, sigma_sec_eff_max = 0;

	min = 0;
	max = 82;
	trouve = -1;

	while (min <= max && trouve == -1)
	{
		if (sigma_ionisation_e_born[(max + min) / 2][0] == E)
		{
			trouve = (max + min) / 2;
			break;
		}
		else if (sigma_ionisation_e_born[(max + min) / 2][0]>E)
			max = (max + min) / 2 - 1;
		else
			min = (max + min) / 2 + 1;
	}

	if (trouve >= 0)
	{
		if (E >= energies_ionisation[4])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][6];
			temp_indice = 4;
		}
		else if (E >= energies_ionisation[3])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3] + sigma_ionisation_e_born[trouve][4];
			temp_indice = 3;
		}
		else if (E >= energies_ionisation[2])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2] + sigma_ionisation_e_born[trouve][3];
			temp_indice = 2;
		}
		else if (E >= energies_ionisation[1])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1] + sigma_ionisation_e_born[trouve][2];
			temp_indice = 1;
		}
		else if (E >= energies_ionisation[0])
		{
			sigma_sec_eff = sigma_ionisation_e_born[trouve][1];
			temp_indice = 0;
		}
		else
		{
			cout << "Erreur : section efficace = 0 GetIndiceEnergieIonisation " << endl;
			exit(1);
		}

		sigma_rand = frand(0, sigma_sec_eff);
		for (i = 1, somme = 0; i <= temp_indice; i++)
		{
			somme += sigma_ionisation_e_born[trouve][i];

			if (somme >= sigma_rand)
				break;
		}
	}
	else
	{
		if (E >= energies_ionisation[4])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][6];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][6];
			temp_indice = 4;
		}
		else if (E >= energies_ionisation[3])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3] + sigma_ionisation_e_born[min][4];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3] + sigma_ionisation_e_born[max][4];
			temp_indice = 3;
		}
		else if (E >= energies_ionisation[2])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2] + sigma_ionisation_e_born[min][3];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2] + sigma_ionisation_e_born[max][3];
			temp_indice = 2;
		}
		else if (E >= energies_ionisation[1])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1] + sigma_ionisation_e_born[min][2];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1] + sigma_ionisation_e_born[max][2];
			temp_indice = 1;
		}
		else if (E >= energies_ionisation[0])
		{
			sigma_sec_eff_min = sigma_ionisation_e_born[min][1];
			sigma_sec_eff_max = sigma_ionisation_e_born[max][1];
			temp_indice = 0;
		}
		else
		{
			cout << "Erreur : section efficace = 0  GetIndiceEnergieIonisation  " << endl;
			exit(1);
		}

		sigma_rand = frand(0, Interpolation(sigma_ionisation_e_born[max][0], sigma_sec_eff_max, sigma_ionisation_e_born[min][0], sigma_sec_eff_min, E));

		for (i = 1, somme = 0; i <= temp_indice; i++)
		{
			somme += Interpolation(sigma_ionisation_e_born[max][0], sigma_ionisation_e_born[max][i], sigma_ionisation_e_born[min][0], sigma_ionisation_e_born[min][i], E);

			if (somme >= sigma_rand)
				break;
		}


	}


	return (i - 1);

}

double QuadInterpolation(double e11, double e12, double e21, double e22, double xs11, double xs12, double xs21, double xs22, double t1, double t2, double t, double e)
{
	double interpolatedvalue1 = (xs11 == xs12) ? xs11 : Interpolation(e11, xs11, e12, xs12, e);
	double interpolatedvalue2 = (xs21 == xs22) ? xs21 : Interpolation(e21, xs21, e22, xs22, e);
	double value = (interpolatedvalue1 == interpolatedvalue2) ? interpolatedvalue1 : Interpolation(t1, interpolatedvalue1, t2, interpolatedvalue2, t);

	return value;
}

void GetDirectionElectronSecondaire(const double & Ti, const double & Ts, double & cosTheta, double & phi)
{
	phi = 2. * PI * UniformRand();

	if (Ts < 50.)
		cosTheta = (2.*UniformRand()) - 1.;
	else if (Ts <= 200.)
	{
		if (UniformRand() <= 0.1) 
			cosTheta = (2.*UniformRand()) - 1.;
		else 
			cosTheta = UniformRand()*(sqrt(2.) / 2.);
	}
	else
	{
		double sin2O = (1. - Ts / Ti) / (1. + Ts / (2.*electron_mass_c2));
		cosTheta = sqrt(1. - sin2O);
	}
}

double GetEnergieElectronSecondaire(double Ti, int indice)
{
	double maximumEnergyTransfer = 0.;

	if ((Ti + energies_ionisation[indice]) / 2. > Ti)
		maximumEnergyTransfer = Ti;
	else
		maximumEnergyTransfer = (Ti + energies_ionisation[indice]) / 2.;

	double crossSectionMaximum = 0.;

	double minEnergy = energies_ionisation[indice];
	double maxEnergy = maximumEnergyTransfer;
	int nEnergySteps = 50;

	double value = minEnergy;
	double stpEnergy(std::pow(maxEnergy / value, 1. / static_cast<double>(nEnergySteps - 1)));
	int step = nEnergySteps;

	double differentialCrossSection;


	int count_nb = 0;

	while (step>0 && count_nb <5)
	{
		step--;
		
		differentialCrossSection = DifferentialCrossSection(Ti, value, indice);

		if (differentialCrossSection >= crossSectionMaximum)
			crossSectionMaximum = differentialCrossSection;
		else
			count_nb++;

		value *= stpEnergy;

	}

	double secondaryElectronKineticEnergy = 0.;

	do
	{
		secondaryElectronKineticEnergy = frand(1e-5, 1.) * (maximumEnergyTransfer - energies_ionisation[indice]);

	} while (UniformRand()*crossSectionMaximum > DifferentialCrossSection(Ti, (secondaryElectronKineticEnergy + energies_ionisation[indice]), indice));

	
	return secondaryElectronKineticEnergy;

}

void GetIndiceSigmaDiff(double Ti, int & min_inf, int & min_sup, int & max_inf, int & max_sup)
{
	int min = 0, max = 82, trouve = -1;

	while (min <= max && trouve == -1)
	{
		if (Tab_Idx[(min + max) / 2][0] == Ti)
		{
			trouve = (min + max) / 2;
			break;
		}
		else if (Tab_Idx[(min + max) / 2][0] < Ti)
			min = ((min + max) / 2) + 1;
		else
			max = ((min + max) / 2) - 1;
	}


	if (trouve >= 0)
	{
		min_inf = (int)Tab_Idx[trouve][1];
		min_sup = (int)Tab_Idx[trouve][2];
		max_inf = (int)Tab_Idx[trouve][1];
		max_sup = (int)Tab_Idx[trouve][2];
	}
	else
	{
		min_inf = (int)Tab_Idx[max][1];
		min_sup = (int)Tab_Idx[max][2];
		max_inf = (int)Tab_Idx[min][1];
		max_sup = (int)Tab_Idx[min][2];
	}

	

}

double DifferentialCrossSection(double Ti, double Tt, int Indice)
{
	double Tt11, Tt12, Tt21, Tt22, xs11, xs12, xs21, xs22, Ti1, Ti2;
	int min_inf, min_sup, max_inf, max_sup, trouve, min, max;

	//Si Energie transfere < plus petite energie eau quitter
	if (Tt  < energies_ionisation[Indice])
		return 0;

	//Chercher encadrement de Ti (2 valeurs 2 lignes) Ti1, Ti2;
	GetIndiceSigmaDiff(Ti, min_inf, min_sup, max_inf, max_sup);

	if (min_inf == max_inf)
	{
		Ti1 = Ti2 = sigmadiff_e_born[min_inf][0]; //valeur exacte de Ti qui est dans sigmadiff_e_born
	}
	else
	{
		Ti1 = sigmadiff_e_born[min_sup][0];
		Ti2 = sigmadiff_e_born[max_inf][0];
	}

	

	//Chercher encadrement de Tt (4 valeurs 2 col pour chaque ligne) Tt11, Tt12, Tt21, Tt22
	min = min_inf;
	max = min_sup;
	trouve = -1;



	while (min <= max && trouve == -1)
	{
		if (Tt == sigmadiff_e_born[(min + max) / 2][1])
		{
			trouve = (min + max) / 2;
			break;
		}
		else if (Tt > sigmadiff_e_born[(min + max) / 2][1])
			min = (min + max) / 2 + 1;
		else
			max = (min + max) / 2 - 1;
	}

	if (min == max || trouve >= 0)
	{
		Tt11 = Tt12 = sigmadiff_e_born[trouve][1];
		xs11 = xs12 = sigmadiff_e_born[trouve][Indice + 2];
	}
	else
	{
		Tt11 = sigmadiff_e_born[max][1];
		Tt12 = sigmadiff_e_born[min][1];
		xs11 = sigmadiff_e_born[max][Indice + 2];
		xs12 = sigmadiff_e_born[min][Indice + 2];
	}

	min = max_inf;
	max = max_sup;
	trouve = -1;

	while (min <= max && trouve == -1)
	{
		if (Tt == sigmadiff_e_born[(min + max) / 2][1])
		{
			trouve = (min + max) / 2;
			break;
		}
		else if (Tt > sigmadiff_e_born[(min + max) / 2][1])
			min = (min + max) / 2 + 1;
		else
			max = (min + max) / 2 - 1;
	}

	if (min == max || trouve >= 0)
	{
		Tt21 = Tt22 = sigmadiff_e_born[trouve][1];
		xs21 = xs22 = sigmadiff_e_born[trouve][Indice + 2];
	}
	else
	{
		Tt21 = sigmadiff_e_born[max][1];
		Tt22 = sigmadiff_e_born[min][1];
		xs21 = sigmadiff_e_born[max][Indice + 2];
		xs22 = sigmadiff_e_born[min][Indice + 2];
	}

	double section_efficace = QuadInterpolation(Tt11, Tt12, Tt21, Tt22, xs11, xs12, xs21, xs22, Ti1, Ti2, Ti, Tt);

	

	return section_efficace;
}
